function vars = updatePosteriorXZ(mix, data, vars)


K = mix.ncentres;
[N D] = size(data.Y);

m = mix.varposterior.m;
W = mix.varposterior.W;
v = mix.varposterior.v;

invC = data.invC;
invCy = data.invCy;

vW = zeros(D, D, K);
vWmt = zeros(K, D);
for k = 1:K
    vW(:, :, k) = v(k)*W(:, :, k);
    vWmt(k, :) = vW(:, :, k)*m(k, :)';
end

E = zeros(D, D, N, K);
R = zeros(N, D, K);
RRt = zeros(D, D, N, K);
logdetE = zeros(N, K);

if strcmp(mix.covar_type, 'full') || strcmp(data.uncertainty_type, 'full')
    % compute E, R, logdetE
    for n = 1:N
        for k = 1:K
            E(:, :, n, k) = inv(invC(:, :, n) + vW(:, :, k));
            R(n, :, k) = E(:, :, n, k)*(invCy(n, :)' + vWmt(k, :)');

            logdetE(n, k) = logdet_chol(E(:, :, n, k));
        end
    end
else
    % faster than repmat
    indsN = false(D, D, N);
    indsNK = false(D, D, N, K);
    indsK = false(D, D, K);
    for d = 1:D
        indsN(d, d, :) = true;
        indsNK(d, d, :, :) = true;
        indsK(d, d, :) = true;
    end
    
    % diagonal update E
    invdiagC = reshape(invC(indsN), D, N);
    diagE = zeros(D, N, K);
    for k = 1:K
        diagE(:, :, k) = 1 ./ (invdiagC + diag(vW(:, :, k))*ones(1, N));
    end
    E(indsNK) = diagE;
    
    % diagonal update R
    for k = 1:K
        R(:, :, k) = (invCy + ones(N, 1)*vWmt(k, :)).*diagE(:, :, k)';
    end

    % diagonal update logdet(E)
    logdetE = reshape(sum(log(diagE), 1), N, K);
end

% compute RRt
for k = 1:K
    RRt(:, :, :, k) = computeCrossProducts(R(:, :, k));
end

vars.E = E;
vars.R = R;
vars.RRt = RRt;
vars.logdetE = logdetE;
